First, we load the data into the workspace.
Following is a summary of the dataset.
The analysis in the frequency domain reveals a strong harmonics with a period of 24 hours. Therefore, it is reasonable to
freq.plot(checkin.NY,title="New York City")
The distribution of the check-in categories accross one period (24 hours).
time.distribution.plot(checkin.NY,title="New York City")
The distribution of the check-in categories in 24 hours in radial plots.
time.radial.plot(checkin.NY)
have a look at a intuitive spatiotemporal visualizaiton of the data…
basemap = map.plot("../../global/data/shapefiles",
"NYC_borough_boundaries_WGS84",
alpha=0.1,size=0.3,color="grey")
saveGIF({
point.animation.plot(checkin.NY,title="New York City",
basemap=basemap,axis.size=18,
title.size=24,legend.size=18,plot.size=3)
}, interval = 0.5, movie.name = "spatiotemporal.gif",
ani.width = 1500, ani.height = 1200)
What will the statistics say when we neglect the temporal factors? First of all, if we try to find spatial clusters based on differnt checkin categories, the distribution of the founded clusters are quite differnt. It indicates that each category has differnt correlations with the geographic space.
if(!exists("basemap")){
basemap <- map.plot("../../global/data/shapefiles",
"NYC_borough_boundaries_WGS84",
alpha=0.1,size=0.3,color="grey")
}
## OGR data source with driver: ESRI Shapefile
## Source: "../../global/data/shapefiles", layer: "NYC_borough_boundaries_WGS84"
## with 5 features and 4 fields
## Feature type: wkbMultiPolygon with 2 dimensions
plots <- lapply(split(checkin.NY,checkin.NY$cate_l1),function(ci){
# find out clusters for the type of category
clusters = spatial.clustering(ci)
centers = clusters[["centers"]]
points = clusters[["point.unique"]]
wss = clusters[["wss"]]
pct = clusters[["pct"]]
lbl = paste(nrow(centers),"clusters;\nWSS:",
formatC(wss,digits=2,format="f"),
"(",format.percent(pct),")")
# add basic points
gg.map <- point.plot(points,x="lon.x",y="lat.x",alpha=0.05, basemap=basemap)
# add cluster information
gg.map <- point.plot(centers, x="lon.center",y="lat.center",
color= "cid.ordered", color.concrete=FALSE,
point.size = log(centers$size,5),alpha = 0.7,
basemap=gg.map,
xlim=range(checkin.NY$lon,finite=TRUE),
ylim=range(checkin.NY$lat,finite=TRUE))
# some plot configuration
gg.map <- gg.map + ggtitle(ci[1,"cate_l1"]) +
theme(legend.position="none") +
geom_text(aes(x = -74.13, y = 40.87), label = lbl, size=2)
})
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
# png("categorized_clustering.png",width=15*ppi,height=6*ppi,res=ppi)
grid.arrange(plots[[1]],plots[[2]],plots[[3]],plots[[4]],
plots[[5]],plots[[6]],plots[[7]],plots[[8]],
plots[[9]],plots[[10]],nrow=2, ncol=5)
# dev.off()
We could also see the most domimant categories in each part of the entire city.
# prepare polygon data
SPDF = readOGR(dsn = "../../global/data/shapefiles", layer = "NYC_zipcode")
## OGR data source with driver: ESRI Shapefile
## Source: "../../global/data/shapefiles", layer: "NYC_zipcode"
## with 236 features and 6 fields
## Feature type: wkbPolygon with 2 dimensions
# intersect the checkin data with the each postal code region
checkin.NY = point.in.poly(checkin.NY, SPDF, copy.attr="POSTAL")
## [1] "Spatial data has been prepared."
## [1] "The information from the polygon has been assigned to the points."
# get the distribution of categories in each postal code region
cate.distr=cate.distr.in.poly(checkin.NY)
## [1] "The distribution of category has been assigned to each polygon."
# and use that information to implement the SPDF
SPDF = poly.with.category(SPDF, cate.distr=cate.distr, poly.attr="POSTAL")
# plot
mapdf=df.from.spdf(SPDF)
# mapdf$density=apply(mapdf,1,function(i){i[i["cate.dom"]]})
# mapdf$density=as.numeric(formatC(mapdf$density,digits=1,format = "f"))
png("main_category.png",height=6*ppi, width=16*ppi,res=ppi)
grid.arrange(
map.plot(mapdf=mapdf,
more.aes=aes_string(fill="Category.1st"),
color="grey",size=0.3,alpha=0.9)+
theme_bw()+
xlab("")+ylab("")+
theme(legend.title = element_text(size=8),
legend.text = element_text(size = 8)),
map.plot(mapdf=mapdf,
more.aes=aes_string(fill="Category.2nd"),
color="grey",size=0.3,alpha=0.7)+
theme_bw()+
xlab("")+ylab("")+
theme(legend.title = element_text(size=8),
legend.text = element_text(size = 8)),
nrow=1, ncol=2, widths=c(1,1) )
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
dev.off()
## pdf
## 2
cate.conds = xtabs(~conds+cate_l1, data=checkin.NY)
#prop.table(cate.conds, 1) # row percentages
#prop.table(cate.conds, 2) # column percentages
fit <- ca(cate.conds)
#print(fit) # basic results
summary(fit) # extended results
##
## Principal inertias (eigenvalues):
##
## dim value % cum% scree plot
## 1 0.003244 52.2 52.2 *************************
## 2 0.002082 33.5 85.8 ****************
## 3 0.000511 8.2 94.0 ****
## 4 0.000147 2.4 96.4 *
## 5 7.6e-050 1.2 97.6 *
## 6 6.7e-050 1.1 98.7
## 7 5.1e-050 0.8 99.5
## 8 2.4e-050 0.4 99.9
## 9 8e-06000 0.1 100.0
## -------- -----
## Total: 0.006210 100.0
##
##
## Rows:
## name mass qlt inr k=1 cor ctr k=2 cor ctr
## 1 | Cler | 501 779 88 | 22 452 76 | -19 327 86 |
## 2 | Fog | 2 634 31 | -265 628 37 | -26 6 1 |
## 3 | Haze | 22 975 178 | -200 798 272 | 94 177 94 |
## 4 | HvyR | 5 917 75 | 208 421 61 | 226 496 111 |
## 5 | HvyS | 0 116 5 | -170 20 0 | 377 96 1 |
## 6 | LgFR | 0 714 29 | -590 712 40 | 28 2 0 |
## 7 | LghR | 33 947 88 | 77 362 61 | 98 586 154 |
## 8 | LghS | 7 246 21 | -62 203 8 | 29 44 3 |
## 9 | MstC | 80 738 38 | -31 319 23 | -35 420 47 |
## 10 | Ovrc | 232 957 107 | -25 217 45 | 46 740 237 |
## 11 | PrtC | 55 592 41 | 38 315 25 | -36 277 34 |
## 12 | Rain | 10 943 67 | 163 632 81 | 114 312 62 |
## 13 | SctC | 49 874 97 | -59 287 53 | -85 587 169 |
## 14 | Snow | 4 850 134 | -421 849 218 | 4 0 0 |
##
## Columns:
## name mass qlt inr k=1 cor ctr k=2 cor ctr
## 1 | ArtE | 105 763 76 | 49 526 76 | 33 237 54 |
## 2 | CllU | 17 677 27 | -75 580 30 | 31 97 8 |
## 3 | Evnt | 5 557 20 | -24 23 1 | -113 534 32 |
## 4 | Food | 294 688 48 | 18 301 28 | -20 388 56 |
## 5 | NghS | 122 955 206 | 93 826 325 | 37 129 80 |
## 6 | OtdR | 69 890 200 | -9 4 2 | -127 886 528 |
## 7 | PrOP | 98 986 272 | -123 878 457 | 43 108 87 |
## 8 | Rsdn | 27 352 42 | 23 54 4 | 54 298 38 |
## 9 | ShpS | 156 661 49 | -29 430 41 | -21 230 34 |
## 10 | TrvT | 107 791 59 | -33 318 36 | 40 473 84 |
#plot(fit) # symmetric map
plot(fit, mass = TRUE, contrib = "absolute", map ="rowgreen",
arrows = c(TRUE, FALSE)) # asymmetric map
Under the assumption \(H=i\) is independent from \(W=j\) \[P(C=k|H=i,W=j)=\frac{P(C=k,H=i,W=j)}{P(H=i,W=j)}=\frac{P(H=i,W=j|C=k)*P(C=k)}{P(H=i)*P(W=j)} (1) \]
since \(H=i\) is independent from \(W=j\), \[Exp[P(H=i,W=j|C=k)]=P(H=i|C=k)*P(W=j|C=k) (2)\]
therefore, \[Exp[P(C=k|H=i,W=j)]=Exp[ \frac{P(H=i,W=j|C=k)*P(C=k)} {P(H=i)*P(W=j)}] \\\ =\frac{P(H=i|C=k)*P(W=j|C=k)*P(C=k)}{P(H=i)*P(W=j)} \\\ =\frac{\frac{P(H=i,C=k)}{P(C=k)}*\frac{P(W=j,C=k)}{P(C=k)}*P(C=k)}{P(H=i)*P(W=j)} \\\ =\frac{P(C=k|H=i)*P(H=i)*P(C=k|W=j)*P(W=j)}{P(H=i)*P(W=j)*P(C=k)} \\\ =\frac{P(C=k|H=i)*P(C=k|W=j)}{P(C=k)} (3)\]
\[P_{u}(C=k|H=i)=\frac{\Phi_{u}(C=k,H=i)}{\Phi_{u} (H=i)} (4)\]
get.temporal.impact <- function(dataframe,hour){
# dataframe.in.hour = checkin.single[which(checkin.single$hour==hour),]
dataframe.in.hour = dataframe[which(dataframe$hour==hour),]
phi.h = nrow(dataframe.in.hour)
list.category = split(dataframe.in.hour, dataframe.in.hour$cate_l2)
sapply(list.category, function(i){
nrow(i)/phi.h
})
}
\[P_{u}(C=k|W=j)=\frac{Intercept(C=k,W=j)}{\sum Intercept(C,W=j)} (5)\]
get.meteorological.impact <- function(fit,conds){
conds.id = which(fit[["rownames"]]==conds)
ref.vec = fit[["rowcoord"]][conds.id,1:8]
cate.all = fit[["colcoord"]][,1:8]
intercepts = apply(cate.all, 1, function(x){
(x[1]*ref.vec[1] + x[2]*ref.vec[2] + x[3]*ref.vec[3] + x[4]*ref.vec[4] +
x[5]*ref.vec[5] + x[6]*ref.vec[6] + x[7]*ref.vec[7] + x[8]*ref.vec[8] ) /
(ref.vec[1]^2 + ref.vec[2]^2 + ref.vec[3]^2 + ref.vec[4]^2 +
ref.vec[5]^2 + ref.vec[6]^2 + ref.vec[7]^2 + ref.vec[8]^2 )
} )
# vec = intercepts / sum(intercepts) !!!wrong!!! intercetp can be negative
vec = (intercepts - min(intercepts)) / (max(intercepts)-min(intercepts)) # scale into [0,1]
names(vec) = fit[["colnames"]]
vec
}
get.meteorological.impact2 <- function(dataframe,conds){
dataframe.in.conds = dataframe[which(dataframe$conds==conds),]
phi.c = nrow(dataframe.in.conds)
list.category = split(dataframe.in.conds, dataframe.in.conds$cate_l2)
sapply(list.category, function(i){
nrow(i)/phi.c
})
}
\[P_{u}^{*} (C=k|W=j)= w_{j}*[P_{u}(C=k|W=j)-\bar P_{u}]+\bar P_{u}\]
get.weather.weight <- function(fit){
conds.all = fit[["rowcoord"]][,1:2]
mag = apply(conds.all, 1, function(x){
sqrt( (x[1]^2+x[2]^2) )
} )
mag / sum(mag)
}
get.weighted.meteorological.impact <- function(fit, conds, weight){
# cate.conds = xtabs(~conds+cate_l2, data=dataframe)
# fit <- ca(cate.conds)
unweighted = get.meteorological.impact(fit,conds)
# weights = get.weather.weight(fit)
# conds.id = which(fit[["rownames"]]==conds)
# vec = weights[conds.id] * (unweighted-mean(unweighted)) + mean(unweighted)
vec = weight * (unweighted-mean(unweighted)) + mean(unweighted)
names(vec) = fit[["colnames"]]
vec
}
get.weighted.meteorological.impact2 <- function(dataframe, conds, weight){
unweighted = get.meteorological.impact2(dataframe, conds)
weight * (unweighted-mean(unweighted)) + mean(unweighted)
}
\[P(C=k)=\frac{\Phi_{u} (C=k) }{\Phi_{u} }\]
get.denominator <- function(dataframe){
phi.h = nrow(dataframe)
list.category = split(dataframe, dataframe$cate_l2)
denominator = sapply(list.category, function(i){
nrow(i)/phi.h
})
denominator
}
\[E[P(C=k|H=i,W=j)]=\frac{P(C=k|H=i)*P(C=k|W=j)}{P(C=k)}\]
update @2014.09.04
p.k has a relatively high performance, multipling it with others do not bringer a higher one; instead, dividing does.update @2014.09.05
get.meteorological.impact(), I used intercept/sum(intercept), where intercept can be negative, and the sum(intercept) can be nonsense. Now it is fixed to (intercepts - min(intercepts)) / (max(intercepts)-min(intercepts)) to scale into [0,1].p.kj is calculated by just using the probability distribution, then the equation p.kij = p.ki * p.kj / p.k is correct (which confirmes our derivation);p.kj is calculated by personal CA, it should be noticed that the scaled interception does not describe the probability of C=k under the condition W=j. Instead, it describes the negative/positive impacts. Hence, the conditional probability should be the scaled interception multipled by the priori p.k. Therefore, here the final result should be p.kij = p.ki * p.kj.p.j, the above two ways have very similar results (which is expected!).